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ABSTRACT 

Non-thermal pressure in galaxy clusters leads to underestimation of the mass of galaxy clus¬ 
ters based on hydrostatic equilibrium with thermal gas pressure. This occurs even for dy¬ 
namically relaxed clusters that are used for calibrating the mass-observable scaling relations. 
We show that the analytical model for non-thermal pressure developed in |Shi & Komatsu] 
(20141 can correct for this so-called ‘hydrostatic mass bias’, if most of the non-thermal pres¬ 
sure comes from bulk and turbulent motions of gas in the intracluster medium. Our correction 
works for the sample average irrespective of the mass estimation method, or the dynamical 
state of the clusters. This makes it possible to correct for the bias in the hydrostatic mass es¬ 
timates from X-ray surface brightness and the Sunyaev-Zel’dovich observations that will be 
available for clusters in a wide range of redshifts and dynamical states. 


Key words: galaxies; clusters: general - galaxies: clusters: intracluster medium - cosmology: 
observations - methods; analytical - methods; numerical 


1 INTRODUCTION 


Galaxy clusters promise great statistical power as a low red- 
shift cosmological probe, if their masses are estimated accurately 
( [Allen, Evrard & Mantz|20fT| l. A bias in the mass estimation poses 
the most serious challenge to the success of cluster cosmology. For 
example, the recently reported tension in cosmological inferences 
from the cosmic microwave background (CMB) and galaxy clus¬ 
ters measured by the Planck satellite, the so-called ‘Planck CMB- 
Cluster tension’ ( [Planck Collaboration V.|20l3] [Planck Collabora-j 
[tion XXIV.[2015[ (, may point to new physics, but it may also reflect 
the urgent need for better understanding and correction for system¬ 
atic biases (see e.g. von^^erLindMe£^J20M Leistedt, Peiris &[ 


[Verde[2014[[Donahue et al.[2014[[Israel et al.[26l5[l.' 

The most prominent astrophysical bias has to do with the 
assumption of hydrostatic equilibrium (HSE) of the intracluster 
medium (ICM) inside the gravitational potential. HSE no longer 
holds in the presence of non-thermal sources of pressure. Hydro- 
dynamical simulations show significant pressure associated with 
bulk and turbulent motions of the gas, which leads to a ‘HSE mass 
bias’ of order ten percent or more in the mass estimates of indi- 


vidual clusters (e.g., Kay et al.|2004|[Rasia et al.|2006 
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As the amplitude of non-thermal pressure can hardly be mea¬ 
sured directly by observations, one has to resort to theoretical es¬ 
timations to correct for this bias. To this end, we have developed 
a physically motivated analytical model for non-thermal pressure 
in the ICM ( [Shi & Komatsu|20T4l SK14 hereafter). Given the ac¬ 
cretion history of the cluster, the model predicts the amplitude of 
non-thermal pressure and its radial, mass, and redshift dependen¬ 
cies. The model successfully reproduces the average non-thermal 
pressure profiles seen in cosmological hydrodynamical simulations 
( [Battaglia et al.[[?012 i, the average thermal pressure profile mea¬ 
sured by Planck lAmaud et al.|2010| (, as well as non-thermal pres¬ 
sure profiles of individual clusters in hydrodynamical simulations 
( [Shi et al.[2015l SKNN15 hereafter). 


In this paper, we study the ability of the SK14 model in cor¬ 
recting the HSE mass estimation by testing it on simulated clusters 
in a set of high-resolution cosmological hydrodynamical simula¬ 
tions ( [Nelson, Lau & Nagai|2014l l, with special attention paid to 
the dynamical state of galaxy clusters. 


The dynamical state of a cluster influences its mass estimation 
in two ways. Eirst, the more dynamically relaxed a cluster is, the 
smaller its non-thermal pressure (e.g. SK14; [Nelson et ^[2014[ l, 
and thus the smaller its HSE mass bias. Second, the cluster profiles 
are smoother when a cluster is dynamically more relaxed, which 
gives a better precision in the mass estimation. 


In the early days, HSE mass estimates were applied to clus¬ 
ters irrespective of their dynamical state (e.g. [Reiprich & Bohringer[ 
\2002\. These studies use X-ray surface brightness data and make 
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the assumption that the ICM is isothermal. After the recognition 
that gas temperature decreases systematically with increasing ra¬ 
dius outside the cluster core region, precise HSE mass estimates 
require also spatially-resolved, high signal-to-noise gas tempera¬ 
ture measurements, which are harder to obtain especially out to 
high redshifts. 

The current cluster mass estimates for a large sample of clus¬ 
ters are mostly performed statistically (rather than individually) 
using scaling relations between the mass and spatially-averaged 
observables, such as the mean X-ray temperature, luminosity, or 


their combination, and the Sunyaev-Zel’dovich (SZ) effect ( |Sun-| 
lyaev & Zeldovich|1970||1972^ averaged over certain cluster radii 


(e.g. 

Vikhlinin et al.||2009 |Andersson et al.||2011 Bender et al. 

|20T4 

IPlanck Collaboration XX.|2014| see also|Giodini et al.|2013 


for a review). In this context, HSE mass estimates are usually per¬ 
formed only on a small sample of dynamically relaxed clusters 
where the bias from non-thermal pressure is minimized, and are 
used for observationally calibrating the scaling relations with the 
hope that the adopted scaling relation is robust against varying dy¬ 
namical states (which holds true for some scaling relations, see 
|Kravtsov, Vikhlinin & Nagai|2006^ . 

However, even for relaxed clusters, non-thermal pressure is 
expected to be non-negligible because non-thermal motions are 
continuously sourced by the mass accretion and are never fully 
thermalized, especially at the cluster outskirts (SK14;|Nelson et al.| 
|2014[ l. Therefore, it is important to study the correction of HSE 
mass bias induced by non-thermal pressure for dynamically relaxed 
clusters. This is one of the primary goals of this paper. 

In recent years, the advances in SZ observations make it pos¬ 
sible to estimate HSE masses using X-ray surface brightness com¬ 
bined with the gas pressure profiles probed by the SZ effect, with¬ 
out the expensive X-ray temperature measurements. This will make 
individual cluster mass estimates for a large number of clusters pos¬ 
sible out to high redshifts, and this time with high precision to large 
radii thanks to the relative insensitivity of the SZ effect to gas den¬ 
sity. In light of these observational advances, another goal of this 
paper is to study the possibility of correcting HSE mass bias for a 
population of clusters in a wide range of redshifts and dynamical 
states. 

The rest of the paper is organized as follows: we introduce our 
method of correcting non-thermal pressure in cluster mass estima¬ 
tion in Sect. where we discover the necessity of regulating the 
input cluster profiles by fitting or smoothing. Then we present four 
methods for regularizing the profiles in Sect. |2.2| and |2.3| and apply 
them to our simulated clusters. After introducing the way we select 
subsample of clusters based on their dynamical state in Sect. |2.4| 
we present our results in Sect. |2.5| and conclude in Sect.[^ In the 
two appendices we study the effect of velocity anisotropy (Ap¬ 
pendix]^ and how an error on the pressure profile propagates to 
an error on the mass (Appendix [B)i. 


2 HSE MASS ESTIMATE AND ITS CORRECTION 
2.1 Method 


The assumption of HSE relates the thermal pressure Pd, that is 
observable from both X-ray and SZ observations to the hydrostatic 
mass Mhse as 


fHHSE(< f) = - 


d dPth 
Gpgadr) dr 


( 1 ) 



Figure 1. Spherically-averaged overdensity profiles (equation]^ with re¬ 
spect to the critical density, A(r), of a simulated cluster. The horizontal dot¬ 
ted lines show A = 500 and 200. The true overdensity profile (black solid 
line) is calculated from the mass profile in the simulation, whereas the es¬ 
timated profiles are calculated from spherically-averaged gas density and 
pressure profiles in the simulation using equations ^ and ^ (red and blue 
lines, respectively). Structures and noises in the gas density and pressure 
profiles are greatly amplified by the derivative and division procedures re¬ 
quired in computing the overdensity profile (dashed lines), suggesting the 
need for fitting (solid lines) or smoothing. 


where Pgas is the gas density. Since the total pressure in the ICM, 
Plot = Fill + Pnth where Fnth is non-thermal pressure support, is 
greater than Eth by definition, Mhse is an underestimate of the true 
mass. We thus define the ‘corrected mass’ Moon 


Mcorr('^ r) — ' 


dPio 


Gpgas(r) dr 


( 2 ) 


This Mcorr is expected to be close to an unbiased estimate of the 
total mass when the non-thermal pressure support is the dominant 
source of the HSE mass bias, as is the case found in hydrodynam- 
ical simulations at least for regions within rjoo (e.g.|Lau, Nagai &| 
|Nelson|2013l|Nelson et al.|2014| l. 

Cluster masses are usually defined with respect to a radius 
within which the mean density equals e.g. A = 500 times the criti¬ 
cal density of the Universe pcut. 


4;rr^ 

M^=-^PonlA. ( 3 ) 

Solving for M^ using equation 0 or 0 is equivalent to find¬ 
ing the corresponding radius from a spherically-averaged ‘over¬ 
density’ profile defined by 


1 3GM(< r) _ 3 1 dP 

Pent 47Tr^ 4;rpcrit rpgodr) dr ’ 


( 4 ) 


with the implicit relation Afr^) = A. Here, P is either Pth or Ptot. 
An error on the pressure gradient due to neglecting Pnth will lead to 
an error on M^ through an incorrect estimation of r^. 

We apply the HSE and corrected mass estimates to a mass- 
limited sample of 65 galaxy clusters with M 500 > 2.2x 10'‘'/!“* Mq 
at redshift z = 0 and their most massive progenitors at z = 0.6 in 


* Unlike for Pth, random motions producing P„th need not be isotropic. 
While we have implicitly assumed the random motions to be isotropic here. 


we investigate the consequence of velocity anisotropy in Appendix [A| 
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the OmegaSOO simulatiorp] ( [Nelson et al.|[2Q14 i, an Eulerian hy¬ 
drodynamic cosmological simulation. Radiative cooling, star for¬ 
mation and feedback are not included in the current runs of the 
OmegaSOO simulation. These additional physics are crucial for re¬ 
producing the observed behaviors of cluster cores. They also af¬ 
fect the detailed shapes of the thermo-dynamical profiles in the 
bulk of ICM, especially for low mass systems. For low mass clus¬ 
ters and galaxy groups, testing non-thermal pressure and cluster 
mass estimations would be better performed with simulations with 
cooling and feedback. However, these additional physics do not 
change the overall shape of the ICM profiles outside the cluster 
core (e.g. |Nagai, Kravtsov & Vikhlinin|2007| l. In particular, the ve¬ 
locity dispersion of the gas and fraction of non-thermal pressure 
at around rsoo are consistent to typically <10% between simula¬ 
tions with and without cooling and feedback ( jNagai et al.||2013l 
[Nelson, Lau & Nagai||2014) even for clusters with relatively low 
masses M 200 ~ 2x lO*'* Mq (corresponding to M 500 ~ 10*'* Mq) 
( [Battaglia et al.|20T^ . Therefore, non-radiative simulations are ad¬ 
equate for our purpose. We refer to Sect. 2 of SKNN15 and ref¬ 
erences therein for further descriptions of the simulation and the 
cluster sample used in this paper. Small clumps have been removed 
when constructing the one-dimensional cluster profiles by exclud¬ 
ing high-density tail in the probability distribution of gas densities 
using the [Zhuravleva et ^ ( [2013| l method. 

We first compute the overdensity profile A(r) from simulated 
profiles with equation (j^, taking F = F^h for Mhse and P = Pth/(1“ 
/nth) for Afcorr- We then use equation 0 to compute Msoo- Here, /„th 
is the ‘non-thermal fraction’ defined by /„th = Pnth/CAh + ^nth)- We 
compute /nth of each cluster using the SK14 model and the mass 
accretion history of the cluster measured from the simulation (see 
SKNN15 for details). 

Numerically calculating a cluster mass using equation 0 is 
not trivial, as taking derivative of pressure and dividing by density 
would amplify small-scale structures and noise. If we apply no reg¬ 
ulation of the simulated profiles other than a cubic-spline interpo¬ 
lation when computing the derivative of the pressure profile, the re¬ 
sulting overdensity profile will be far from smooth (dashed lines in 
Fig. despite that the gas density and pressure profiles are rather 
smooth after clumps are removed. In the previous X-ray studies, 
the observed gas density and temperature profiles are often fit to 
parametric models to reduce this noise (e.g. Pointecouteau, Arnaud| 


|& Pratt|2005||Vikhlinin et al.|2006[[Schmidt & Allen|2007| l. Non- 

parametric smoothing methods have also been adopted in simula¬ 
tion studies jLau, Kravtsov & Nagai|200^|Nelson et al.|2014^ . 

Here we use and compare four approaches: one based on cubic 
spline smoothing, and the others based on fitting pressure and den¬ 
sity to dilferent parametrized models. For the parametrized models, 
we further compare two fitting ranges. 


2.2 Cubic spline smoothing 

First, we compute the overdensity profile (equation [^ from 
simulated gas density and pressure profiles with cubic spline inter¬ 
polation applied to pressure when taking a derivative, dFjdr. We 
then smooth the computed log A(log r) profile with a cubic-spline 
smoothing algorithm described in |Dierclgt]jl993^. The algorithm 


^ The sample is chosen with a mass limit of M 2 oom > 6 X M© at 

z = 0, with M 2 oom being the mass enclosed in a radius within which the 
average density is 200 times the mean matter density of the universe. 


fits the data points with a piece-wise cubic spline with an increas¬ 
ing number of nodes, until the input constraint on the value of x~ 
statistic is satisfied. We use the constraint that divided by the 
number of degrees of freedom is less than 0.04. This constraint is 
chosen to minimize the variance of MnsE/iHcorr while not biasing 
the sample mean. 

When estimating M 500 with the smoothing method, a large ra¬ 
dial range of data around the true r^go is needed to guarantee that 
the estimated A(r) profile covers A = 500. Therefore, we apply this 
method only to rather relaxed clusters using data on the radial range 
of [0.1, 1.5] 


2.3 Parametrized models 

Next, we fit parametrized models to the pressure and gas den¬ 
sity profiles and estimate masses. We fit pressure rather than tem¬ 
perature because pressure is the quantity that directly enters the 
HSE and corrected mass estimation. Also, thermal pressure is the 
quantity that the thermal Sunyaev-ZeTdovich (SZ) effect ^Sunyaev[ 
|& Zeldovich|197()[|1972^ directly probes, so fitting gas density and 
pressure may be a better option when applying to combined SZ and 
X-ray data. 

The data points used in the fit are equally spaced logarithmic 
bins with uniform weighting on the radial range of either [0.1, 1] 
'^500 [^'1’ '^ 500 - These are the typical radial ranges of the 

currently available data on well-resolved clusters used for the mass 
calibration. 


2.3.1 Model ‘VW 


We first use the model presented in |Vikhlinin et ah] ([2006|l, 
which is based on an isothermal (3 model ( [Cavaliere & Fuscc^ 
|Femiano|1978^ with modifications capturing the cool cores and the 
observed steepening of the gas density profile at large radii. Here 
we ignore the modification for the cool cores, both because cool 
cores are non-existent in our simulated clusters, and that our focus 
is on the outer region of the clusters, especially around rsoo where 
we determine their masses. The resulting simplified model is: 


Pgas(f) Of 


1 1 
(1 -Hr2/r2)T (1 -I-ri'Vr")^ ’ 


(5) 


and 


Pir) Of Pgas(r) 


(rlr,r 
(1 +P’//f)f/'’ ■ 


( 6 ) 


Following jVikhlinin et al.H2006| l, we fix = 3. After this, this 
model still has 10 parameters (including 2 for normalization), and is 
highly nonlinear, suggesting the possible existence of local minima 
in the function. 

For all the parametrized models, we fit the simulated gas den¬ 
sity and pressure profiles to the models by searching for the global 
minimum of the x^ function in a wide parameter range. Here, we 
set 0 < e < 5, 0 < yS < 10,-1 < a < 5, 0 < h < 10, 0 < c < 10, and 
and r, being less than 1 Mpc/h. 


2.3.2 Model ‘VG’ 

The Vikhlinin model was designed to fit gas density and tem¬ 
perature of the ICM. It is then not surprising that the fitting form for 
the pressure (equation [^ is rather redundant. Here, instead, we use 


© 0000 RAS, MNRAS 000,(3-?? 





































4 Xun Shi, Eiichiro Komatsu, Daisuke Nagai and Erwin T. Lau 


a generalized Navarro-Frenk-White (GNFW) profile i Zhao|1996) 
for pressure, as suggested by|Nagai, Kravtsov & Vikhlinin 1 2007 1 : 


P{r) oc 


r’' [l + (r/f-p)F“j 


aifi-y) 


(V) 


This simpler form already has enough flexibility to describe the 
shape of the observed thermal pressure profiles (e.g. |Arnaud et al.| 
|2010| >, as well as our simulated and Ptot- Flere we further fix 
7 = 0 which provides an adequate fit to all our simulated clusters. 
For the GNFW parameters, we search for the best fit values within 
0.5 <a<2, and Vp between 10 kpc/h and 1 Mpc/h. 

For the gas density we still use the Vikhlinin formula (equation 
[^. The result is a 9-parameter model. 


2.3.3 Model ‘NG’ 


An ideal parametrized model fits the data with enough flex¬ 
ibility with a minimum number of fitting parameters. The way to 
achieve this is to utilize existing knowledge or well-founded as¬ 
sumptions on the underlying physical system that the data repre¬ 
sent. 

Flere we try to do so by assuming that the total mass density 
profile follows an NFW profile ( [Navarro, Frenk & White|[l996| ). 
The gravitational acceleration, g(r), is then given by 


^ ^ GM« r) (1 -H r/r,) ln(l + r/r,) - r/r, 
gir) = -5- oc — 


rHl + r/r,) 


( 8 ) 


where is the scale radius of NFW profile, whose value we search 
for between 10 kpc/h and 1 Mpc/h. 

For P,h and Ptot we continue to use GNFW (equation]^. The 
gas density then follows from P{r) and g(r) as 


Pgas(r) — 


1 dP(r) 
g(r) dr 


(1 -l-r/rj)[7-l-/?(r/rp)F“j 

[(1 -I- r/rj)ln(l -I- r/r,) - r/r,] r)'-‘ [l -I- (r/rp)'/o']“‘^ 


(9) 


After setting 7 = 0, the model has 6 parameters. 

The overdensity profiles constructed with this model for one 
representative cluster in the sample are shown as the colored solid 
lines in Fig.[T] For this cluster, the fits to the gas density and total 
pressure profiles (i.e., corrected for non-thermal pressure) give an 
overdensity profile (blue solid line) that agrees very well with the 
true one (black solid line). This suggests that at least for this cluster, 
the NFW assumption works and this model provides a good fit to 
the simulated profile. 


2.4 Dynamical state 

Sourced by the growth of cluster via mergers and accretion, 
non-thermal pressure correlates strongly with the mass growth his¬ 
tory. The latter also correlates with the dynamical state of the clus¬ 
ter, but the degree of correlation depends on the particular classifi¬ 
cation of dynamical states. Here, we try to mimic the observational 
procedure and determine the dynamical state of the clusters based 
on their X-ray morphologies. This allows us to estimate the level of 
non-thermal pressure for the ‘relaxed’ clusters that are used to cal¬ 
ibrate the mass-observable scaling relations. We compare our clas¬ 
sification to that using an automatic method developed by |Mantz| 
jet al.| ( |20r5t in Appendix]^ 

We use mock Chandra X-ray images in a single projec¬ 
tion, which are produced using the pipeline developed in [Nagai,] 



Figure 2. Mock X-ray images of four clusters representative of different 
dynamical states: ‘very relaxed’ (upper left, cluster number 1 in Fig.|^, 
‘rather relaxed’ (upper right, cluster number 5 in Fig.|^, ‘less disturbed’ 
(lower left), and ‘disturbed’ (lower right). The white dashed circles show 
the positions of rsoo. The exposure time used for the mocks is 100 kilo- 
seconds. 


[Vikhlinin & Rravtsov[j200'7| (. In particular, we examine (1) whether 
the X-ray image shows a single, distinguished cluster core with lit¬ 
tle displacement with respect to the bulk of the ICM; (2) whether 
the cluster core region appears relaxed with round or elliptical con¬ 
tours; (3) whether there are substructures or clear disturbances in 
the ICM between cluster core region and rsoo; and (4) whether there 
are substructures or clear disturbances between [1, 1.5] rsoo that 
will influence the pressure gradient estimation at rjoo- 

Only 1 out of the 65 clusters appears very relaxed according 
to all the above criteria. If we relax the criteria to allow for the ex¬ 
istence of small concentrated substructures that can be removed by 
the clump-removing method, then 3 more clusters can be classified 
as ‘very relaxed’. We conshuct larger sub-samples of clusters by 
further relaxing the criteria: a ‘rather relaxed’ sample of 14 clusters 
which have no or only slight disturbances in the core region and the 
global ICM, and contain no or only small concentrated substruc¬ 
tures within and around r;oo', and a ‘less disturbed’ sample of 32 
clusters (including the ‘rather relaxed’ clusters) which have single, 
distinguished cluster core and no significant substructure between 
cluster core region and r^oo- Fig. shows mock Chandra X-ray im¬ 
ages of representative clusters of these samples. 


2.5 Results 

2.5.7 Influence of mass estimation methods on the HSE mass 
bias 

In the upper panel of Fig. we show the HSE mass biases of 
the 14 ‘rather relaxed’ clusters. Given the same input gas density 
and pressure profiles, the HSE masses estimated using different ap¬ 
proaches (Sect. and [2.3| ( vary at the level of 5-10% for all 14 
clusters. None of the four approaches systematically give signifi¬ 
cantly smaller or larger HSE mass biases, suggesting that they are 
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Table 1. Numerical values of Bias ± Scatter in mass estimations with the fitting method ‘NG’ shown by the red symbols and error bars in Fig.[^ The values 
for the progenitors at z = 0.6 are also shown in the last column. 



fitting range 

rather relaxed (14/65) 

less disturbed (32/65) 

mass-limited (65/65) 

progenitors at z=0.6 (65/65) 

ln(MHSE/Wtiiie) 

[0.1,l]r<- 
[0.1, 1.5] 

-0.09 ±0.14 
-0.16 ±0.08 

-0.14 ±0.23 
-0.21 ±0.14 

-0.23 ±0.31 
-0.24 ±0.15 

-0.24 ± 0.26 

-0.26 ±0.15 

ln( Mqqxi / ^true ) 

[0.1. 1] r™ 
[0.1, 1.5] 

-0.01 ±0.18 

-0.07 ±0.11 

-0.02 ± 0.27 
-0.04 ±0.17 

-0.04 ± 0.32 
-0.03 ± 0.20 

-0.07 ± 0.28 

-0.05 ±0.18 
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Figure 3. HSE mass biases at r 5 oo (upper panel) and their correction with 
the simulated non-thermal pressure (middle panel) and the modeled non- 
thermal pressure (lower panel) for the ‘rather relaxed’ clusters. The 14 
clusters are ordered by their degree of relaxation, ranging from very re¬ 
laxed on the left (number 1) to slightly disturbed on the right (number 14). 
The masses are estimated using profiles between 0.1 and 1.5 with four 
methods mentioned in the text: a spline smoothing method (S) and three 
parametrized model fitting methods (NG, VV and VG). The averaged mass 
biases for the sample before and after correction are given in Fig. and 
Table 1. 


comparable in a statistical sense. Even for the few most relaxed 
clusters (cluster number 1-4), there exist 5-10% HSE mass biases, 
confirming the need for non-thermal pressure correction even for 
the relaxed cluster sample typically used for calibrating the scaling 
relation. 

Fitting to different radial ranges also introduces 5-10% dif¬ 
ferences in the HSE mass estimates (Fig.|^. Fitting to larger radii 
more likely gives larger HSE mass biases. This is because the in¬ 
crease of non-thermal pressure with radius causes additional steep¬ 
ening of Fth compared to Ftot, so when data from larger radii are 



lnMHSE[fit over 0.1-1 r-*So']-luM,rue 

Figure 4. Effect of fitting range on the HSE mass bias. HSE mass biases 
at r 5 oo of the ‘rather relaxed’ 14 clusters estimated using the pai'ametrized 
fitting method ‘NG’ are shown. 



NG VV VG S NG VV VG NG VV VG 


Figure 5. Averaged HSE (open markers) and corrected (filled markers) 
mass biases at rsoo and their scatter for cluster samples of different relax¬ 
ation states: top 14 (‘rather relaxed’) and 32 (‘less disturbed’), and the full 
mass limited 65 clusters (‘all’). The fitting / smoothing method is indicated 
by the label on the x-axis, with symbols identical to that in Fig.[^ Fitting 
methods are performed using profiles on two different radial ranges: [0.1, 1 ] 
^500 (^ii'^l^s) and [0.1, 1.5] (squares). The numerical values are given 
in Table 1. 


considered in the fit, larger HSE mass biases are favored given 
the limited flexibility of any fitting formula. Being a weak effect, 
this systematic trend is visible only for relatively relaxed clusters 

(Fig.|g. 
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2.5.2 Correction for the HSE mass bias 

How well can we correct for the HSE mass bias with the 
SK14 analytical model for non-thermal pressure? The lower panel 
of Fig. shows the corrected masses for the ‘rather relaxed’ clus¬ 
ters. The correction for individual clusters is not perfect, which 
is also the case if we correct for the masses using the simulated 
Enth (middle panel). The imperfection of HSE mass correction with 
simulated Enth is also noticed in |Nelson et al.|p014| ). This imper¬ 
fection is most likely due to asphericity and the small scale struc¬ 
ture of the ICM as well as residual accelerations which also breaks 
the HSE assumption ( |Lau, Nagai & Nelson||20l3) ). However, the 
average HSE mass bias for the 14 clusters is greatly reduced af¬ 
ter correction, both with simulated and modeled irrespective 
of which smoothing/fitting method is used. This is more apparent 
from Fig.|^ where we show the averaged HSE (open markers) and 
corrected (filled markers) mass biases and their scatter. Remark¬ 
ably, for the ‘less disturbed’ and even the full mass-limited sample, 
both at z=0 (5th column of Table 1) and a high redshift of z=0.6 
(6th column), the biases are successfully reduced to be well be¬ 
low 10% irrespective of the smoothing/fitting method or the fitting 
range (Fig.|^ Table 1). This suggests that, the HSE mass biases 
related to non-thermal pressure can be corrected robustly for the 
sample average with the SK14 model for clusters in a wide range 
of dynamical states and redshifts, which will be obtained with com¬ 
bined X-ray surface brightness and SZ observations. 


3 CONCLUSIONS AND DISCUSSIONS 

There is no universal value for the HSE mass bias due to non- 
thermal pressure. Not only does it depend on masses, radii, and 
redshifts (SK14), but it also depends on the particular fitting or 
smoothing method used to compute the HSE mass from the ob¬ 
served profiles, as well as on the radial range of the data used for 
the fitting. Thus, to evaluate how much non-thermal pressure con¬ 
tributes to the ‘Planck CMB-cluster tension’, the particular method 
used by the Planck team needs to be tested e.g. on simulated cluster 
samples. 

Non-thermal pressure biases the masses, hence the mass- 
observable scaling relations, even for the most relaxed clusters in 
the mass-limited sample in the simulation. This poses a challenge to 
the commonly adopted approach of calibrating the scaling relation 
on a sample of relaxed clusters and applying it to the full cluster 
sample at all dynamical states. Correction for non-thermal pressure 
is needed for the calibration of the scaling relation on relaxed clus¬ 
ters too because non-thermal motions are continuously sourced by 
the mass accretion and are never fully thermalized, especially at the 
cluster outskirts. 

The SK14 model is able to correct the non-thermal pressure 
biases in the mass-limited sample and its dynamically more relaxed 
subsamples, irrespective of the dynamical state, the smoothing or 
the fitting method used for the mass estimation, or the radial range 
of data used in the fitting. Aided by this correction, the HSE mass 
estimate may no longer be limited by non-thermal pressure or dy¬ 
namical state of the cluster. This opens up a possibility for more 
accurate cluster mass determinations using combined X-ray sur¬ 
face brightness and SZ observations. Compared to cluster masses 
estimated using weak gravitational lensing, the bias-corrected HSE 
masses are more precise, i.e. have smaller statistical uncertainties. 
Nevertheless, it will be important to check the consistency of both 
methods. 


Gas motions in the ICM will be probed by the upcoming 


ASTRO-H mission (jTakahashi et al.l2012b through the shifts and 

broadening of the X-ray spectral lines 

ZuHone, Markevitch &| 

|Zhuravleva|2015||Ota, Nagai & Lau|2015 

1 . However, the measure- 


ments will be limited to a small number of clusters and small radii 
(r ^ r 25 oo) where the random motions are mild. Although the mea¬ 
surements cannot be directly used to correct mass estimations at 
rsoo, they will test our current understanding on the nature of gas 
motions in the ICM. 

The necessary input of the SK14 model is the mass ac¬ 
cretion history of clusters. Although not directly observable, the 
mass accretion history is known in a sample-averaged sense from 
dark matter-only cosmological simulations (e.g. |Zhao et al.|2009| >, 
though its cosmology dependence awaits more detailed simula¬ 
tion studies. Our simple correction method reduces the mean of 
the mass biases but not yet the scatter, which means no advantage 
will be gained by correcting HSE mass biases for individual clus¬ 
ters at this stage. If the residual acceleration, asphericity and small 
scale structures of the ICM, which are the most likely sources of the 
scatter, are better understood in the future, then knowledge of non- 
thermal pressure in individual clusters will be very helpful in get¬ 
ting accurate individual masses. The mass accretion history and the 
non-thermal pressure for individual clusters may also be obtainable 
by exploiting the connection between the mass density profile and 
the mass accretion history ( [Ludlow et al.|2013||Diemer & Kravtsov| 
|2014| (. Namely, first compute a mass profile from a cluster without 
correcting for non-thermal pressure, and then compute the mass 
accretion history from the mass profile, compute non-thermal pres¬ 
sure, and recompute the mass profile. Then iterate until the mass 
profile converges. We leave this for future work. 
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APPENDIX A: VELOCITY ANISOTROPY 

In the absence of streaming motions, the steady-state radial 
Jeans equation gives (see e.g. |Lau, Nagai & Nelson|2013^ 

GM{< r) _ \ d (Ah + PgasO-?) (2o-2 - cr^) 

^ - 7 » (Al) 

^ Pgas or r 



Figure Al. Difference in the Mcorr obtained by assuming and not assuming 
j3 = 0 when correcting for non-thermal pressure is smaller than the correc¬ 
tion itself, Mtrue ” ^HSEi for the majority of clusters in the simulation. The 
color scale shows p: red for bigger p and blue for smaller p. The stars show 
the subsample of relaxed clusters, while the circles show the rest. Only 10 
out of 65 clusters show |Mcorr,^ - Mcorr^=ol > A^tme “ A^HSE- 


where crj is the radial non-thermal velocity dispersion, and cr^ = 
is the tangential non-thermal velocity dispersion. 

Let us define the velocity anisotropy parameter, p, as 


O't 


(A2) 


This par amet er vanishes for isotropic velocity distribution, and 
equation jAl| reduces to equation l|^ with Rn,i, = PgasO'^th- Here, 
= (07 +o- 2)/3 is the averaged one-dimensional non-thermal 
velocity dispersion squared. 

We calculate the influence of a non-zero anisotropy parame¬ 
ter on the HSE mass estimate. Inserting equation ( |A2| l into equa¬ 
tion l|Al[l, we obtain 


M(< r) = - 


fl In Pgas 
31nr 




Cpgas dr 

_ 

dlnr 


3-2/J G 
ain(3-2y3) 


fllnr 


+ 2li 


(A3) 


When ft t 0, the difference between the mass estimates without 
assuming ft = 0, Mp, and assuming ft = 0, is given by 


-I- 


'■) - ^/^=o(< r)] = 

2ft aino-g„ 3 a In (3-2/1) 
3-2ft dlnr 3-2ft dlnr 


d In Pgas 
d\nr 


+ 


3-2ft 


(A4) 


The first two terms on the r.h.s. dominate; however, they cancel 
when Pgas has a logarithmic slope of -3, which happens to coincide 
with the outer slope of an NEW profile. At rsoo, Pgas typically has a 
slightly shallower logarithmic slope than -3, but still letting these 
terms cancel to a large extent. Therefore, the difference in mass cor¬ 
rections with and without including velocity anisotropy is expected 
to be smaller than the mass correction itself, which is 


G 31n p„as 

-^ [Mcorr« r) - MhsE« C] = 

ro-2„ dlnr 


+ 


dlnr 


(A5) 


In this Appendix, we use the cubic spline smoothing method 
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/nth 


Figure Bl. Correlation between non-thermal pressure and the HSE mass 
bias at rsoo- The open and filled symbols show MnsE/j^tme from the ‘NG’ 
parametric fitting and the cubic spline smoothing methods, respectively, as 
a function of the non-thermal fraction, /nth, at rsoo- The stars show the sub¬ 
sample of relaxed clusters, while the triangles show the rest. The dashed 
and solid lines show the slopes of -1 and -2.5 to guide the eye. 


to regularize the simulated cluster profiles in the radial range of 
shows that the difference between assum- 


A1 


[0.01, 3] r'™. Fig. _ 

ing and not assuming /I = 0 is indeed smaller than the correction 
itself, i.e., - Mhse. for the majority of clusters. Even in the 

cases where the difference is significant, the mass estimated by tak¬ 
ing into account the velocity anisotropy is not always closer to the 
true mass. Thus including velocity anisotropy does not give signifi¬ 
cantly improved results. This, together with that j8(r) is not observ¬ 
able, makes us conclude that the simple mass correction given in 
equation 0 would already suffice. 


APPENDIX B: FROM NON-THERMAL PRESSURE TO 
MASS BIAS 

In this Appendix, we show how an error on the pressure pro¬ 
file propagates to an error on the mass. According to equation 0. 
an error on the pressure gradient will lead to an error on A(r) as 
din A = -6ln(dP/dr). Since the error on the estimation of 
is related to the error on the A(r) profile through the chain rule 
dln(r) = dlnAfr) ((91nA/31nr)“*, the cluster mass will be miss- 
estimated by 


dln(MA) = 3dln(rA) =. (Bl) 

The slope of A(r) is approximately -2 (±0.2) at a radius of r^oo 
for an NEW mass density profile with a concentration parameter 
typical of clusters. Hence, (51n(il/A) ~ 36\n(dP/dr) /2. 

The relation between the change in the pressure gradient 
d In (dP/dr) and the non-thermal fraction /„th is: 


din 




dr j 


= ln(l -/nth) + ln 


1 ± 


/ <91nP,o, y 

\ dlnr / 


ain(l-/„,h) 

51nr 


(B2) 


method used in regularizing the input cluster profiles. We show in 
Fig. |B1 [ that |dln(MA)| estimated with both the ‘NG’ fitting and the 
cubic spline smoothing methods correlate positively with /„*, and 
the fitting method yields a significantly larger HSE mass bias since 
we fit out to a large radius of here. 

If the HSE mass estimation is performed at the true ta (as e.g. 
in Sect. 6 of SK14), the mass bias is dln(MA) = Sln(dP/dr), i.e., 
about 2/3 of the above value. 


APPENDIX C: COMPARISON OF DYNAMICAL STATE 
CLASSIFICATIONS 


Recently, |Mantz et al.| ( |20l5| l developed an automatic method 
for classifying the relaxation state of galaxy clusters. Their classi¬ 
fication is based on the symmetry (‘s’), peakiness (‘p’), and align¬ 
ment (‘a’) of their X-ray morphologies, similar to the factors we 
consider for our classification (Sect. [2!4^ . Here we compare the 
two classification methods on our simulated mass-limited cluster 
sample at z = 0 using mock X-ray images with 100 kilo-seconds 
exposure time (Fig. |Cl^ . 

In |Mantz et al.| l 2015^ classification, a cluster is categorized as 
relaxed when a s-p-a criterion of ^ > 0.87, p > -0.82, and a > 1.00 
(dashed lines in Fig. |Cl^ is satisfied in more than 50% of the cases 
of its bootstrap analysis. This results in a relaxed cluster fraction of 
about 16%, slightly smaller than the fraction (21%) of our ‘rather 
relaxed’ clusters (blue circles and the green star in|Mantz et al.| 
( |2015^ ). In general, more relaxed clusters in our simulated sam¬ 
ple according to our classification also have higher s-p-a values, 
showing overall agreement of the two classification methods. The 
specific distribution of s-p-a values of our simulated cluster sam¬ 
ple is slightly shifted compared to the distribution of morphologi¬ 
cal values of the Chandra sample (Fig. 8, |Mantz et al.| ( [20T5] l): our 
simulated clusters have systematically higher alignment values, and 
lower peakiness values. This is likely due to the inability of non- 
radiative simulations in reproducing the observed properties of the 
cluster core. After taking the corresponding shifted s-p-a criterion 
into consideration, the distribution of s-p-a values of our ‘rather 
relaxed’ clusters aligns rather well with that of the relaxed clusters 
as classified with|Mantz et al.|(|2015|) method. 


In the limit that /„th is small and constant with radius, the rela¬ 
tion is linear: (51n(AfA) ~ -3/nth/2. In practice, is radial depen¬ 
dent (|Shi & Komatsu|2014|>, and this relation also depends on the 
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Figure Cl. Distributions of |Mantz et al.|(2015^ morphological values for 
the simulated mass-limited cluster sample at z = 0. The s-p-a cuts defining 
the relaxed sub-sample given in |Mantz et al.|(2Q15^ are shown as dashed 
lines. In comparison, galaxy clusters classified into different dynamical state 
categories using method given in this paper are indicated by different mark¬ 
ers: ‘very relaxed’ (green star), ‘rather relaxed’ (blue circles and the green 
star), ‘less disturbed’ (magenta triangles, blue circles and the green star), 
and ‘disturbed’ (red crosses). 
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